RNA-Seq-Based WGCNA and Association Analysis Reveal the Key Regulatory Module and Genes Responding to Salt Stress in Wheat Roots

Soil salinization is the main abiotic stressor faced by crops. An improved understanding of the transcriptional response to salt stress in roots, the organ directly exposed to a high salinity environment, can inform breeding strategies to enhance tolerance and increase crop yield. Here, RNA-sequencing was performed on the roots of salt-tolerant wheat breeding line CH7034 at 0, 1, 6, 24, and 48 h after NaCl treatment. Based on transcriptome data, a weighted gene co-expression network analysis (WGCNA) was constructed, and five gene co-expression modules were obtained, of which the blue module was correlated with the time course of salt stress at 1 and 48 h. Two GO terms containing 249 differentially expressed genes (DEGs) related to osmotic stress response and salt-stress response were enriched in the blue module. These DEGs were subsequently used for association analysis with a set of wheat germplasm resources, and the results showed that four genes, namely a Walls Are Thin 1-related gene (TaWAT), an aquaporin gene (TaAQP), a glutathione S-transfer gene (TaGST), and a zinc finger gene (TaZFP), were associated with the root salt-tolerance phenotype. Using the four candidate genes as hub genes, a co-expression network was constructed with another 20 DEGs with edge weights greater than 0.6. The network showed that TaWAT and TaAQP were mainly co-expressed with fifteen interacting DEGs 1 h after salt treatment, while TaGST and TaZFP were mainly co-expressed with five interacting DEGs 48 h after salt treatment. This study provides key modules and candidate genes for understanding the salt-stress response mechanism in wheat roots.


Introduction
Soil salinization is a universal environmental issue that results in detrimental abiotic stress in agriculture, which inhibits crop growth and reduces grain yield [1]. Approximately 1 billion hectares of land worldwide are currently salinized, representing about 7% of the Earth's land surface and 33% of irrigated lands [2].With the continuous growth of the global population, it is necessary to fully utilize salinized land to increase crop yields and meet the growing demand for food.As the most widely planted crop in the world, common wheat (Triticum aestivum L., 2n = 6x = 42; AABBDD) contributes about one-fifth of the total calories and proteins of daily dietary intake [3].To ensure food security, it is crucial to explore salt-tolerance genes from germplasm and to utilize them to enhance the adaptability of wheat varieties in saline soil.
For common wheat, which has an enormous and complex hexaploid genome, transcriptome sequencing is a rapid and effective method for mining candidate genes [4][5][6], especially in studies involving a stress response.Based on RNA-sequencing (RNA-Seq) data, many important biological or abiotic stress response genes have been further identified, such as the leaf rust-resistance gene Lr47 [7] and the cold tolerance gene TaSnRK1α [8].RNA-Seq has also been applied in salt tolerance research.In Kharchia Local, a well-known salt-tolerant wheat germplasm resource in India, 2495 differentially expressed genes (DEGs) in the roots [9] and 3197 DEGs in the flag leaves [10] were identified by comparing the transcriptomes of seedlings in the NaCl-treated group and the control group, and 22 saltinduced unigenes were selected and validated through qRT-PCR.Similar analysis on Saudi salt-tolerant wheat germplasm resource No. 193 showed that 5829 genes were differentially expressed in the roots, whereas 3495 genes were identified in the shoots, of which eight DEGs were validated [11].Moreover, RNA-Seq data from other salt-tolerant germplasm, such as the Chinese cultivars Qingmai 6 [12], Xiaoyan 60 [13], Jimai 22 [14], and Cangmai 6005 [15], as well as Iranian cultivar Arg [16], also revealed many DEGs responding to salt stress, helping to further identify salt-tolerance genes in wheat.However, most of the previously reported RNA-Seq analyses have been based solely on DEGs and use simple qRT-PCR to validate candidate genes, which cannot uncover more useful information.The salt tolerance of wheat is a complex quantitative trait controlled by multiple genes, making it more suitable to use network analysis methods to identify important salt-tolerance response pathways and strongly related genes.
In model plants, such as Arabidopsis, seedlings respond strongly at the molecular level in the early stage of salt stress.As the first organ to be exposed to stress, the roots directly face osmotic stress due to high salt concentrations in the soil [17] and suffer from salt stress, as well as the entrance of Na + through a nonselective cation channel or other means [18].Subsequently, pathways such as signal transduction and transcriptional regulation respond, initiating the plant's salt-tolerance mechanism [19,20].With the development of high-throughput sequencing technology, weighted gene co-expression network analysis (WGCNA) is becoming increasingly mature [21].This method can obtain a series of biologically significant co-expression modules, and has been successfully used to analyze the transcriptome data of multiple samples in rice [22] and maize [23], identifying key regulatory pathways and genes responsive to salt stress.WGCNA was also adopted in a recent study on wheat that reported the expression patterns and the time course of root tissue from salt-tolerant wheat line ST9644 under NaCl stress [24].Here, we performed RNA-Seq on the roots of salt-tolerant wheat breeding line CH7034 [25] after 0, 1, 6, 24, and 48 h of NaCl treatment, and performed WGCNA to screen for important gene co-expression modules related to salt stress.Then, we used salt-tolerant phenotype data from a set of wheat germplasm to perform association analyses on the genes in the selected co-expression modules to identify key salt-tolerance genes.

RNA-Seq Data Evaluation
The transcriptome sequencing of 15 root samples of the salt-tolerant breeding line CH7034 resulted in 166.70 Gb clean reads.Each sample obtained clean data above 9.54 Gb with Q30-based percentages greater than 93.60% and GC percentages ranging from 52.06 to 54.05% (Table S1), indicating that the RNA-Seq data were of high quality for further analysis.Subsequently, 452,237,372 clean reads were mapped on the Chinese Spring reference genome RefSeq v1.0, and the uni-transcripts were annotated as known genes.A total of 56,420 genes with an average fragments per kilobase of exon model per million mapped fragments (FPKM) of >1 in at least one treatment (1, 6, 24, and 48 h) were considered expressed genes.
Differential expression analysis identified 6571, 19,414, 20,273, and 22,105 DEGs at 1, 6, 24, and 48 h of NaCl treatment, respectively, compared with the transcriptome data at 0 h (Figure 1a).The number of DEGs upregulated by salt stress was greater than that of downregulated DEGs (Figure 1a).Overall, 30,369 DEGs were upregulated or downregulated in at least one treatment, and 4017 DEGs appeared at each stress time point (Figure 1b).

Gene Co-Expression Modules Responding to Salt Stress in Roots
WGCNA was performed on the 30,369 DEGs screened above, and 4202 genes were clustered and divided into five modules with different colors based on their expression patterns (Figure 2a).Among them, the gray module contained two genes that could not be assigned to any module, and the remaining four modules, namely turquoise, blue, brown, and yellow, contained 2305, 1706, 118, and 71 genes, respectively (Table S2).The Pearson correlation coefficient (r > 0.5) and the significance of the p value (p < 0.05) were used as the criteria to measure the relationship between module and time course of salt stress.The results revealed a significant correlation between the blue module and time points 0, 1, and 48 h, with "blue module-1 h" showing a positive correlation (r = 0.53; p = 0.04) and "blue module-48 h" showing a negative correlation (r = −0.65;p = 0.009).Moreover, a positive and negative correlation was shown between the yellow module and time points 0 h (r = 0.69; p = 0.005) and 48 h (r = −0.58;p = 0.02), respectively.The brown module was negatively correlated with 0 h (r = −0.76;p = 0.001) and positively correlated with 1 h (r = 0.57; p = 0.03).The turquoise module was negatively correlated with 0 h (r = −0.69;p = 0.005) and positively correlated with 24 h (r = 0.53; p = 0.04), and the gray module was not related to any time point (Figure 2b).
Compared with the other modules, the blue module was associated with more salttreatment time points (0, 1, 48 h) and showed a strong correlation at 48 h.Therefore, all genes in the blue module were used for subsequent Gene ontology (GO) enrichment analysis.

The GO Enrichment of the Blue Module
A total of 455 terms were enriched from the blue module, of which 302 involved biological processes, 86 involved molecular functions, and 67 involved cellular components (Table S3).The 20 terms with the highest −log 10 (FDR, false discovery rate) values, including 11 cellular component terms and nine biological process terms, are shown in Figure 3.The cellular component terms contained Cytosolic ribosome, Nucleolus, and Ribosomal subunit, among others, indicating a close correlation with gene transcription and protein synthesis.Furthermore, the biological process terms included Response to osmotic stress (GO:0006970) and Response to salt stress (GO:0009651), the two key GO terms related to the salt-stress response (Figure 3).Interestingly, all 236 DEGs of the GO:0009651 term were included in the 249 DEGs of the GO:0006970 term; hence, we focused on these 249 DEGs, which could be related to osmotic/salt stress (Table S4).

Gene Co-Expression Networks Responding to Short-Term Salt Stress in Roots
Twenty genes with the highest edge weight co-expressed with the four candidate genes were selected from the 249 DEGs to construct an interaction network.The results showed that 3, 6, 21, and 22 genes interacted with TaGST, TaWAT, TaAQP, and TaZFP, respectively (Figure 6, Table S7).These interacting genes were divided into two categories.One category contained regulatory genes, including TraesCS3A02G413700 encoding Broad Complex, Tramtrack, and Bric-A-Brac (BTB)/POX Virus and Zinc Finger (POZ) protein, and TraesCS6A02G301900 and TraesCS6B02G331000, encoding ethylene-responsive factor (ERF).The other category contained the following structural genes: TraesCS3A02G313600, TraesCS5B02G511800, and TraesCS5B02G529300, encoding kinase; TraesCS3B02G290200 and TraesCS7A02G189500, encoding glycosyltransferase; and TraesCS4B02G304300, encoding a detoxification protein (Table S7).
Among the genes in this network, 17 genes, including TaWAT and TaAQP, were upregulated in the early stage of NaCl treatment, while the remaining seven genes, including TaZFP and TaGST, were upregulated in the later stages of stress, mainly at 48 h (Figure 7).

The Blue Module Is Related to the Time Course of Salt-Stress Response in Wheat Roots
The roots are the first organ exposed to high salt stress.Exploring the key modules and genes in roots is crucial not only for a deeper understanding of salt-tolerance mechanisms but also for salt-tolerance molecular breeding in wheat.Research in model plants has shown that cellular responses can be placed in different phases after salt treatment [26].Early signaling represents changes observed within 5 min of salt stress application in the roots.The transcriptional response from the nucleus is activated, and the growth rate decreases during the stop phase and the quiescent phase within 5 min to 5 h and 5-9 h of salt application, respectively.Later, plants partially recover their growth rate during the growth recovery phase after 9 h of salt application [26].
In this study, we set four time points for the NaCl treatment of roots: 1 h during the stop phase, 6 h during the quiescent phase, and 24 and 48 h during the growth recovery phase.The RNA-Seq results showed that, compared with the untreated control, the number of DEGs increased with the prolongation of treatment time, and there were more upregulated genes than downregulated genes.These DEGs were divided into five modules through WGCNA, of which the blue module was significantly positively correlated with 0 and 1 h and significantly negatively correlated with 48 h, indicating that this module was related to the time course and primarily functioned during the stop and growth recovery phases.Transcriptional changes occurring in the stop phase promote the synthesis of functional proteins, such as biochemical enzymes and ion transporters [27], and activate stress-related ABA or ethylene signaling pathways [19,28].In the later growth recovery phase, genes involved in growth repair begin to be regulated [29].

Four Novel Candidate Genes Responding to Salt Stress Discovered from the Blue Module
Two biological process terms, Response to osmotic stress (GO:0006970) and Response to salt stress (GO:0009651), were enriched in the blue module.Of 249 DEGs, 4, namely TaWAT, TaAQP, TaZFP, and TaGST, from the two GO terms were significantly correlated with salt-tolerance phenotypes in the roots using a set of wheat germplasm.
TaWAT and TaAQP were upregulated in the roots of CH7034 at 1 h of salt treatment.The WAT1 gene encoding the transmembrane protein was involved in auxin transport and homeostasis [30], and was significantly upregulated by salinity stress [31,32].In this study, TaWAT was associated with RSIR-RtL, -RsA, -RT, and -RV, indicating that this gene may ultimately impact the ability of roots to transport auxin and alter growth under salt application, as predicted in another study [33].Aquaporins located on the cell plasma membrane mediate the regulation of root hydraulic conductivity in response to osmotic stress and salt stress [17,34].The overexpression of two wheat AQP genes, TaNIP (GenBank No. DQ530420) [35] and TaAQP8 (GenBank No. HQ650110) [36], significantly enhanced salt tolerance in the transgenic plants of Arabidopsis and tobacco, respectively.The upregulation of TaAQP8 under salt stress involved an ethylene signaling pathway, causing a positive effect and then increased root elongation [36].Moreover, a recent study reported that the AQP activity of wheat varied with the diurnal cycle, which periodically affected the plant's tolerance to salt stress [37].In our case, a novel TaAQP gene was associated with RSIR-RV, and exhibited similar expression patterns 24 and 48 h after salt treatment (Figure 7).
TaZFP and TaGST were mainly upregulated in CH7034 roots 48 h after salt treatment during the growth recovery phase.ZFP, located in the nucleus, is a type of well-studied transcription factor related to salt stress responsiveness [38,39].ZFP overexpression significantly increased the levels of ascorbic acid (AsA) in Arabidopsis and tomato, which strengthened the AsA-mediated reactive oxygen species (ROS)-scavenging capacity and enhanced the salt tolerance of the plants [40].In salt-tolerant wheat cultivar Shanrong No. 3, the ZFP gene TaCHP (GenBank No. GQ379226) was expressed in the roots of seedlings at the three-leaf stage, and its transcript was localized in the cells of the root tip cortex and meristem [41].The TaZFP reported in this study is a novel gene associated with RSIR-RV, and the other gene, TaGST, is related to RSIR-RtL and RSIR-RsA.Overexpressing GST in cotton or soybean improves their salt tolerance [42][43][44].In addition, GST is involved in scavenging ROS induced by salt stress in cells [45].Therefore, TaGST and TaZFP may recover cell growth by alleviating oxidative stress, thereby enhancing the roots' tolerance to salt stress.

A Predicted Gene Co-Expression Network for Salt Tolerance
Based on edge weight values, an interaction network was predicted, setting the four candidate genes as hub genes.After salt stress application to roots, TaWAT and TaAQP were upregulated and promoted the expression of multiple genes involving biochemical reactions.Some of these genes were coding kinase proteins.In the salt-stress response in Arabidopsis, kinase proteins play a crucial role in ion transport and in signal transduction [46,47].Some genes referred to transcription and translation, such as the 60S ribosomal protein L5 (RPL-5), ATP synthase subunit b (ATP5H), and RNA-binding protein (RBP), which are related to salt tolerance in rice or cotton [42][43][44], and some genes referred to the metabolic processes in response to salt stress [48,49], such as the aminotransferase-related protein gene and ketol-acid reductoisomerase NADP(+) (KARI).In addition, two ERF genes, TraesCS6A02G301900 and TraesCS6B02G331000, were also induced.Currently, three salt stress-response ERF genes have been identified in wheat, including two salt-tolerant genes, TaERF1 [50] and TaERF3 [51], as well as one salt-sensitive gene, TaERF4 [52].The two ERFs identified in this study are novel genes related to salt-stress responsiveness.
In the 6 to 48 h after salt stress, TaZFP and TaGST gene activities gradually increased, and multiple genes were induced.TraesCS4B02G304300 encodes a detoxification protein that may alleviate cation toxicity in cells [53,54], and TraesCS2A02G032500 encodes ATP sulfurylase that can chelate SO 4 2− under salt stress [55].TraesCS3A02G413700 encodes a BTB/POZ transcription factor, and TraesCS3B02G290200 and TraesCS7A02G189500 encode glycosyltransferases.Both BTB/POZ and glycosyltransferases have been reported to scavenge ROS, which are induced by salt stress [56][57][58].Therefore, the genes that respond during this period may primarily perform growth recovery functions under salt stress.
However, the blue module was significantly negatively correlated with the 48 h after the salt-stress time course.One possible explanation is that some genes, especially the hub genes TaZFP and/or TaGST, exist as negative regulatory factors in certain salttolerance pathways.Similar results have been reported.Several ZFP genes, such as MtZPT2-2 in Medicago truncatula [59], AZF1 and AZF2 in Arabidopsis [60], and DST [61] and OsZFP185 [62] in rice, can negatively regulate plant salt tolerance by inhibiting the ABA signaling pathway or transport proteins.Moreover, the GST gene AtGSTU17 has also been shown to play a negative role in salt stress tolerance in Arabidopsis [63].Further experiments are needed to confirm this speculation.

Plant Materials and Salt Treatment
Wheat breeding line CH7034 with salt tolerance was used for RNA-Seq and qRT-PCR.A set of wheat germplasm containing 114 varieties (Table S8) was used for salt-tolerance phenotype identification and association analysis.
Sterilized seeds were laid on Petri dishes covered with moist filter paper for 2-3 days.The uniform germinated seeds were then selected and transferred to half-strength Hoagland's culture solution in a growth chamber under a 22/16 • C (day/night) temperature regime and a 16/8 h (light/dark) photoperiod with 60% relative humidity.When the seedlings grew to the three-leaf stage, they were exposed to 250 mM NaCl for salt-stress treatment.

RNA-Seq
The root samples of the CH7034 seedlings were collected with three replications after 0, 1, 6, 24, and 48 h of NaCl stress, then frozen in liquid nitrogen, and stored at −80 • C. Total RNA was extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions, and the cDNA libraries were constructed using TruSeq RNA Sample Preparation Kit v2 (Illumina, San Diego, CA, USA).RNA-Seq was performed on all libraries by Biomarker Technologies Co., Ltd.(Beijing, China) using the HiSeq 4000 platform (Illumina, San Diego, CA, USA).Clean reads were obtained by removing lowquality reads containing adapters or poly-N from the raw data and mapped to the Chinese Spring reference genome RefSeq v1.0 (http://wheat-urgi.ver-sailles.inra.fr/,accessed on 10 March 2023).The transcript level of each gene was measured with FPKM and then normalized by log2.Differential expression analysis between the control (0 h of NaCl treatment) and salt-stress groups (1, 6, 24, and 48 h of NaCl treatment) was performed using DESeq R package [64], with a threshold of |log2FoldChange| ≥ 1 and FDR ≤ 0.01.

Co-Expression Network Construction
Gene co-expression networks were constructed using R (version 4.3.1)software and the WGCNA package (version 1.72) for the RNA-Seq data from the roots.The transcript levels of all of the DEGs were converted into a similarity matrix and transformed to a topological overlap matrix using a parameter β value of 12. Genes with similar expression patterns were categorized into different modules using a bottom-up algorithm with a module minimum size cutoff of 30.
The correlation between module eigengenes and the time points for salt stress was calculated using a Pearson test, and the individual modules with p < 0.05 were considered significantly correlated with the time course.Then, the module with the highest correlation coefficient was used for GO enrichment analysis on the BMKCloud platform (www.biocloud.net,accessed on 24 November 2023), and items with a high correlation with salt tolerance were selected.Genes with a WGCNA edge weight > 0.6 in the selected items were visualized using Cytoscape 3.7.2software.

Association Analysis
RSIR-R was investigated for 114 wheat germplasm resources.Salt stress was imposed for germplasm according to the method described in Section 4.1.A Root Scanner (MicroTek, Shanghai, China) was used to scan the wheat roots after 7 days of H 2 O treatment (CK) and NaCl treatment to obtain the total root length (RtL), total surface area (RsA), total volume (RV), average diameter (RD), root tip number (RT), and root branching number (RF).The relative salt-injury rate of each root phenotype for wheat variety was calculated using the following formula: RSIR-R (%) = (X CK − X NaCl )/X CK × 100%.
Information on sequence variation in the gene coding region was downloaded from the WheatUnion database (http://wheat.cau.edu.cn/WheatUnion/,accessed on 24 November 2023).Using 114 wheat germplasm resources, the correlation between the root-relative saltinjury rate data and the genes screened by WGCNA, the edge weight was evaluated using association analysis in the R (version 4.3.1)software with the GAPIT package (version 3).The correlation was considered significant when −log 10 (p) > 2 (i.e., p < 0.01).

Statistical Analysis
The Origin software (version 3.1) was used to perform the statistical analysis by one-way analysis of variance (ANOVA), and p < 0.05 was considered a statistically significant difference.

Conclusions
Based on RNA-sequencing data performed for the roots of salt-tolerant wheat breeding line CH7034 at 0, 1, 6, 24, and 48 h after NaCl treatment, we conducted WGCNA and obtained five gene co-expression modules, of which the blue module was correlated with the time course of salt stress at 1 and 48 h.Two GO terms containing 249 DEGs related to osmotic stress response and salt-stress response were enriched in the blue module, and 4 DEGs, including TaWAT, TaAQP, TaGST, and TaZFP, were associated with the root salt-tolerance phenotype in a set of wheat germplasm resources.A co-expression network constructed with the four candidate genes and another 20 DEGs showed that TaWAT and TaAQP were mainly co-expressed with 15 interacting DEGs 1 h after salt treatment, while TaGST and TaZFP were mainly co-expressed with 5 interacting DEGs 48 h after salt treatment.

Figure 2 .
Figure 2. Weighted gene co-expression network analysis (WGCNA) in salt-stress response of wheat roots.(a) Cluster dendrogram and module colors of 4202 DEGs.(b) Heatmap for the relationships of modules and time courses.Each cell contains the corresponding correlation and p-value in parentheses.* indicates p < 0.05, and ** indicates p < 0.01.

Figure 3 .
Figure 3. Gene ontology (GO) enrichment analysis of DEGs in the blue-module.The brown and black GO terms belong to the categories of biological process and cellular components, respectively.

Figure 6 .
Figure 6.Co-expression network for four hub genes and 20 DEGs with the highest edge weight values from selected terms in the blue module.The red, green, and blue circles represent hub-genes, regulatory genes, and structural genes, respectively.Arrows indicate the direction of interaction.

Figure 7 .
Figure 7. Heatmap of genes in the co-expression network at 0, 1, 6, 24, and 48 h after salt stress in wheat root.The red, green, and blue genes represent hub-genes, regulatory genes, and structural genes, respectively, and the gene annotations based on Chinese Spring RefSeq v1.1 are listed in parentheses.